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MULTIPLE  BEAMFORMING 
by 

Walter  M.X.  Zimmer 


ABSTRACT 


■  Beamf ormi ng  techniques  are  generally  based  on  the  assumption  that  the 
received  signal  arrives  from  a  single  source.  In  a  realistic  situation, 

however,  there  is  more  than  one  source  present  and  these  sources  can 

interact  with  each  other  in  the  beamformer.  This  paper  describes  a^ 

/-multiple  beamformer  in  which  the  interaction  of  different  sources  is 

included  in  the  algorithm.  The  maximum  likelihood  principle  is  used  to 
find  the  different  source  amplitude  and  bearing  estimates.  The  performance 
of  the  algorithm  is  compared  for  varying  conditions  with  the  Cramer-Rao 
lower  bound.  Two  pre-processors  necessary  and  suitable  to  ensure  globally 
optimized  results  are  derived  from  the  ideas  of  Pisarenko  and  Prony. 
Finally,  a  post-processor  is  proposed  to  evaluate  the  success  of  the 
multiple  beamformer. 


INTRODUCTION 


In  general  beamforming  is  understood  as  a  procedure  for  combining  spatially 
distributed  measurements,  either  in  a  predefined  or  in  an  adaptive  way,  to 
increase  the  sensibility  of  the  receiving  system  in  the  directions  of 
interest. 

The  usual  assumption  made  in  the  implementation  of  a  beamformer  algorithm 
is  that  the  signal  arrives  from  a  single  direction  and  the  interaction  of 
sources  from  other  directions  can  be  neglected.  However,  in  realistic 
situations,  we  have  to  modify  this  assumption  since  we  know  either  that 
there  is  more  than  one  source  or  that  we  have  multipath  arrivals,  therefore 
we  should  include  the  interaction  between  the  sources  or  the  multipath 
arrivals. 

In  this  paper  we  describe  a^method  that  allows  us  to  form  multiple  beams  in 
order  to  steer  the  receiver  toward  two  or  more  sound  sources  simulta¬ 
neously.  The  main  advantage  of  this  procedure  is  that  it  includes  the  full 
interactions  of  the  different  sources  in  the  beamformer.  :  t  > ^  -  s  » 

We  expect  to  gain  the  resolution  of  closely  spaced  sources  without  losing 
too  much  of  the  robust  behaviour  of  the  conventional  beamformer  and  without 
biasing  the  estimate  of  the  quantities  we  want  to  know. 

Statistics  suggest  that  the  Cramer-Rao  lower  bound  (CRLB)  is  a  good  measure 
of  the  robustness  of  a  signal  processor.  If  there  exists  a  method  to 
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achieve  this  bound  then  it  will  be  the  maximum  likelihood  parameter  estima¬ 
tion  technique.  Therefore  it  is  sensible  to  base  the  multiple  beamformer 
on  this  technique  [1  to  5], 

Due  to  the  explicit  modelling  of  reality,  parameter  estimation  techniques 
are  easy  to  understand  and  to  modify.  This  explains  the  widespread  use  of 
such  techniques. 

Why  has  this  rather  simple  approach  not  received  more  attention  in  the 
literature?  The  answer  can  be  found  in  the  fact  that  the  problems  we  are 
dealing  with  generally  do  not  have  unimodal  solutions  and  require  either  an 
exhaustive  search  [5]  or  appropriate  preprocessing.  However,  with  the 
development  of  so-called  high  resolution  methods  suitable  techniques  are 
available  to  achieve  optimized  results. 

The  complete  algorithm  presented  in  this  paper  is  composed  of  four  parts. 
The  essential  components  are  the  multiple  beamformer  that  estimates  the 
amplitudes  of  the  signals  in  closed  form  and  the  fine  bearing  estimator 
necessary  to  calculate  the  correct  phase  relationships  between  the  dif¬ 
ferent  sound  sources.  For  the  fine  bearing  estimator  we  use  an  iterative 
technique  to  find  the  solution.  We  need  a  preprocessor  to  begin  the  itera¬ 
tion.  The  fourth  component  estimates  the  background  noise  level. 

This  paper  first  presents  the  concept  and  algorithm  of  the  multiple  beam- 
former  including  the  fine  bearing  estimator.  Secondly  it  analyses  the 
performance  of  this  kernel,  with  respect  to  the  detection  and  resolution 
capabilities.  The  necessary  preprocessor  will  only  be  sketched  and  some 
possible  implementations  discussed.  The  paper  concludes  with  a  proposal 
for  a  postprocessor. 


1  FORMULATION  OF  THE  BEAMFORMING  PROBLEM 

The  formulation  will  be  restricted  to  passive  sonar  applications,  where  the 
receiver  is  assumed  to  be  a  line  array  with  equidistantly  spaced  hydro¬ 
phones.  Further,  we  assume  that  our  data  are  preprocessed  with  a 
narrowband  filter  bank. 

Therefore  we  have  measurements  with  complex  values  as  a  function  of  space 
and  time: 

Yn.t  »  (i) 


where 

n  is  the  hydrophone  Index  n  =  0,  ...  ,  M-l  and 

t  is  the  sample  index  in  time  t  *  0 . T-l. 


We  wish  to  estimate  the  sound  field  producing  these  hydrophone  measure¬ 
ments  : 


(ak,t»  wk,t) » 


(2) 
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where 

a|(  t  is  the  sound  amplitude, 

wk’t  is  the  source  bearing  angle  and 

k  *  1,  ,  K  is  the  source  index. 

From  the  principle  of  superposition  we  know  that  the  sound  fields  from  dif¬ 
ferent  sources  add  up  linearly  in  space.  Assuming  also  linear  transfer 
characteristics  for  the  hydrophones  and  for  the  narrowband  filter  bank  we 
find  a  linear  relation  between  the  sound  field  and  the  measurements.  This 
means  that  we  can  model  our  measurements  as  follows 

yn,t  =  +n,t(“t)*at  +  wn,t  (3> 


where 


I'n.t^t)  *  (*n,t(“lt) . *n,t(“Kt)) 

is  the  sound  field-measurement  transfer  vector 

at  =  (alt*  •••  »  aKt) 

is  the  soundfield  vector  and 

w„tt  is  the  measurement  error 

*  denotes  the  complex  conjugate  transpose. 

Clearly  there  is  no  unique  solution  to  the  beamforming  problem  as  defined 
by  Eq.  3.  On  the  one  hand  we  have  uncertain  measurements  due  to  noise 
effects;  on  the  other  hand  the  number  of  unknown  parameters  we  wish  to 
estimate  generally  will  not  coincide  with  the  number  of  measurements  we 
have  available. 

Thus,  the  problem  formulation  must  be  restricted  by  introducing  some  sort 
of  a  priori  information  such  as  known  features  of  the  solution,  optimality 
criteria,  or  constraints. 

This  Information  is  usually  introduced  to  reduce  the  beamforming  problem  so 
that  a  useful  solution  can  be  achieved. 


2  DEFINITION  OF  THE  SCENARIO 

It  is  now  necessary  to  define  the  scenario  for  which  the  solution  should  be 
valid.  The  scenario  is  based  on  the  following  assumptions: 

-  The  received  sound  field  consists  of  only  a  small  number  of 
point  sources.  This  excludes  the  presence  of  coloured 
background  noise.  (This  assumption  is  made  for  convenience  but 
In  principle  the  model  can  be  extended  with  terms  describing 
coloured  background  noise.) 

-  The  point  sources  are  In  the  extreme  far  field  so  that  they  can 
be  treated  as  stationary  planewave  signals.  (This  assumption 
is  normally  fullfilled  in  the  case  of  distant  shipping.  The 
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case  of  moving  targets  can  be  included  in  the  model  but  then 
the  computational  workload  will  increase.  Multipath  arrivals 
are  treated  as  different  but  correlated  sources.) 

The  array  speed  is  known.  The  array  speed  is  included  in  the 
algorithm  to  allow  coherent  processing  not  only  in  space  but 
also  in  time. 

The  measurement  errors  and  the  signal  amplitudes  are  random 
variables  which  can  be  approximated  by  gaussian  distributions: 

ak ,t  =  N<ak»  ^k)2  (4) 


wn,t  =  N{0.  o„) 


where 


N(m ,  a2) 


represents  a  complex  valued  gaussian  variable  with  mean  m 

2 

and  variance  o  . 


For  a  line  array  with  equidistantly  spaced  hydrophones  we  can  now  specify 
the  k-the  element  of  the  soundfield-measurement  transfer  vector 


'I'n.tK.t)  58  e 


“k.t* (n+Yt) 


with 


uk.t 


X 


~  c°s(“>k,t) 


(5) 


y 


v-T 


s 


d 


where 

x  is  the  signal  wavelength 
d  is  the  hydrophone  spacing 
v  is  the  array  speed 
Ts  is  the  narrowband  sampling  interval 

Because  of  these  assumptions,  one  must  be  careful  when  applying  tne  results 
to  real-world  situations. 
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3  THE  OPTIMIZATION  CRITERION 

Due  to  measurement  errors  it  is  impossible  to  build  up  a  set  of  consistent 
equations  to  solve  the  beamforming  problem;  but  an  optimal  solution  can  be 
found. 

Which  optimization  criterion  should  we  use?  The  natural  way  to  answer  this 
question  seems  to  be  that  the  optimal  set  of  parameters  6  is  found  when 
the  average  cost  function. 


C(9) 


I  i  C(e- 9)  p(y  ,9 )  dy  de  , 


(6) 


becomes  minimal  [6].  Here  C(6-0)  describes  the  cost  as  a  function  of  the 
estimation  error  9-e  and  ^(y.6)  is  the  joint  probability  density  of 
measurement  y  and  parameter  e. 

Having  no  additional  information  about  the  cost  function  C(8-9)  it  seems 
proper  to  select  constant  cost  for  undesired  deviation  of  the  parameter  6 
from  the  optimal  value  9  thus. 


C(9  -0  )  =  p(  6-8  >  a) 


(7) 


where 


p  (expression)  = 


1  (if  expression  is  true) 


0  (if  expression  is  false) 
is  the  generalized  Kroneker  symbol. 

Minimizing  the  average  cost  yields  the  maximum  a  posteriori  (MAP)  estimate 
defined  by  the  following  equation  [6]: 


—  1o9  p(0  |  y)|  A  =0 

30  10  =  0fi)AP 

Using  the  relation 

p(0,y)  =  p(8  |  y)  p(y )  =  p(y  j  0 )  p(8) 


(8) 


the  MAP  equation  reads 
3 


3Q-  i°g  p(y  |  e)|  A  +l_iogP(e) 

•8  =  6Map  88  8  = 


=  o 


(9) 


10  =  ©map 
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This  equation  shows  how  the  a  priori  density  p(8 )  has  to  be  used  to 
obtain  an  optimal  estimate  of  the  parameter  6.  In  the  case  of  a  constant 
a  priori  density  p(8 )  tne  second  term  disappears  and  the  optimal  estimate 
is  then  called  the  maximum  likelihood  estimate  (MLE) 

=0  (10) 

0  =  0MLE 

In  this  paper  we  assume  that  the  a  priori  density  p(6)  is  constant  so 
that  it  is  sufficient  to  work  only  with  Eq.  10,  simplified  to 

1-1=0,  (ID 

36 


—  '°g  p(y  I  0) 


where 


L  =  log  p(y  |  6)  is  the  log  likelihood  function. 


4  MULTIPLE  BEAMFORMER 

To  define  the  multiple  beamformer  we  must  first  recall  the  basic  relation 
between  the  sound  field  and  the  measurements,  as  given  in  Eq.  3, 

>n,t  =  i'n ,t (“ )"  a  +  "n,t  »  (12) 


where  the  signal  bearing,  and  amplitude  vector,  a^,  are  now  assumed  to 
be  time  independent  (indicated  by  eliminating  the  subscript,  t). 


The  measurement  errors,  wntt,  are  assumed  to  be  independent  gaussian 
variables  in  space  and  time therefore  we  can  write  for  the  log  likelihood 
function,  t,: 


L  =  const  -  2  log  a 


yn,t  -  »ft,t 

a 


where 


l 

n,t 


M-l 

i-  l 

MT  n=0 


T-l 


i 

t=0 


4) 


(13) 


This  symbol,  [  will  be  used  throughout  this  paper. 
n,t 
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With  the  definitions 

b  =  l  yn,t  *n,t 
n  ,t 

*  =  1  ’t’n.t  t'n.t 
n  ,t 

Tr(R)  =  M  •  I  |yn,t|  2 
n,t 

and  after  some  algebraic  manipulations  the  log  likelihood  function,  L,  can 
be  rewritten  as 

L  =  const  -  2  log  a  -  a"2{($a-b)*  $_1(*a-b)  +  Tr(R)-b*  *_1b)  .  (14) 

It  is  obvious  that  the  Eq.  14  is  at  maximum  when 

*a-b  =  0  .  (15) 

This  matrix  equation  can  be  solved  algebraically  thus, 

a  =  *-!b  ,  (16) 

defining  the  multiple  beamformer. 

This  definition  is  well -posed  because  the  matrix  4  includes  all  interac¬ 
tions  between  the  independent  source  directions. 

From  both  Eq.  16  and  the  earlier  definition  of  b  we  can  see  that  this 
multiple  beamformer  processes  the  data  coherently  in  space  and  in  time. 

2 

The  estimation  of  the  variance  a  can  be  found  using  the  equation 

i-  l  I  A  =  0  (17) 

ao  I  o  =  OMLE 


which  yields 


(18) 
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or,  with  Eq.  15, 


o'  =  I  Tr(R)  -  b*  *-ib 


(19) 


Summarizing,  we  can  write  for  the  multiple  beamformer, 
a  =  «_1  b, 

and  for  the  variance  of  the  estimation  error, 

o2  =  ~  Tr(R)  -  b*  $-‘b  .  (20) 

M 

The  residual  log  likelihood  function  is  given  by 

Lr  =  const  -  log  -  {  Tr ( R )  -  b*  «-rb}  .  (21) 

M 


The  multiple  beamformer  as  defined  by  Eqs.  16  or  20  allows  us  to  estimate 
amplitude  and  phase  for  a  finite  number  of  sources  if  the  source  directions 
are  known  in  advance.  In  the  next  section  we  will  weaken  this  condition  by 
requiring  that  the  source  directions  be  known  only  roughly.  Then  a  fine 
bearing  estimator  will  be  used  to  find  the  optimal  bearing  estimates. 


5  TINE  BEARING  ESTIMATOR 

To  estimate  the  fine  bearing  of  the  different  sources  the  residual  log 
likelihood  function,  Lr,  as  given  in  Eq.  21,  can  be  used  as  an  optimization 
measure  since  it  does  not  depend  explicitly  on  the  amplitudes. 

As  usual  we  find  the  desired  bearing  estimates  by  maximizing  this  residual 
log  likelihood,  which  is  also  equivalent  to  the  minimization  of  the 
variance  of  the  amplitude  estimation  error  o2 .  In  this  sense  the  maximum 
likelihood  equation  can  be  written  as 

—  Lrl  =  —  oZ|  =  0  (22) 

3u  lu=UMLE  3u  lu=uMLE 


where  u  is  the  parameter  vector  we  wish  to  estimate,  as  defined  in  Eq.  6: 
2d 

Uk  =  T  - —  COS 
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Unfortunately  Eq.  22  has  no  simple  solution  so  that  we  are  forced  to  use 
some  iterative  techniques. 

Although  more  efficient  techniques  exist,  for  reasons  of  simplicity  we  have 
chosen  the  Newton  procedure  to  solve  Eq.  22. 

As  generally  known,  Eq.  22  is  replaced  in  the  Newton  procedure  by  a  first- 
order  Taylor  approximation  of  the  solution,  which  yields: 


With  the  definitions 


AUj  =  Uj  -  Uj 


(24) 


equation  23  can  be  written  as 

g  =  -  H  Au  ,  (25) 

having  the  solution 

Au  =  -  H-l  g  .  (26) 

Being  close  to  the  MLE  solution  of  Eq.  22,  this  equation  can  be  used  to 
iterate  towards  the  optimal  bearing  estimate. 

We  now  develop  the  expressions  for  the  gradient,  g,  and  the  Hessian,  H. 
Starting  with  Eq.  15 

*a  =  b 


we  obtain 
ab 

3u^ 


and 


k* uj  3Ui/  \3 “  1  / \ 3 U j /  V3Uj/\3Ui/  \3uj  9ui/ 


(27) 


(28) 
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o2  .1 

M 

we  find 

9  i  =  - 

and 


Hij  =  a 


;.  15  and  19 

Tr ( R )  -  b*a 

gL  a  -  b* 

9ui  9u-j 

,*/  32*  \a  -  2/'3a*N)< 

\9Uj  9 uj  '9Uj/ 

(29) 


(30) 


where 


9a  _  9b  9$  g 
Su|  ”  9'uf  ’  3u j 


(31) 


This  completes  the  general  derivation  of  the  terms  necessary  for  the  Newton 
iteration. 

Next  we  have  to  determine  under  which  conditions  the  Newton  iteration  will 
converge  to  the  desired  solution.  The  simplest  way  to  sketch  the  problems 
that  can  occur,  is  to  investigate  the  graphic  representation  of  a  simple 
case. 

Figure  1  shows  the  log  likelihood  function,  the  gradient,  and  the  Hessian 
of  a  single  source  scenario. 

As  indicated,  we  identify  three  areas: 

Area  I  The  global  maximum  can  be  found  by  use  of  some  maximum  gradient 
search  techniques. 

Area  II  The  global  maximum  can  be  found  via  the  Newton  iteration. 

Area  III  The  Hessian  is  negative  definite,  which  is  necessary  for  the 

solution  to  be  a  maximum. 

Within  Area  II  we  see  that  the  iteration  approaches  the  solution  in  an 
alternating  way.  Outside  Area  II  but  within  Area  I  the  Newton  iteration 
diverges  away  from  the  true  solution.  There  is  also  a  singular  behaviour 
of  the  Hessian  at  the  border  of  Area  III. 

The  following  technique  allows  js  to  combine  the  advantages  of  the  Newton 
algorithm  within  Area  II  with  a  simple  gradient  technique.  This  ensures 
convergence  to  the  desired  maximum  inside  Area  I. 
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Considering  the  general  formula  for  the  Newton  algorithm 

Au  =  -  H-i  g  (32) 


we  first  replace  the  regular  inverse  by  the  generalized  inverse,  i.e.. 


H-‘  -  H- 


(33) 


This  is  done  in  order  to  avoid  problems  with  the  singularities  of  the 
Hessian.  Remember  that  the  generalized  inverse  H"  is  any  matrix  satisfying 
the  relation 

H  H-  H  =  H  (34) 


In  the  case  of  a  non-singular  Hessian  the  generalized  inverse  is  equal  to 
the  regular  inverse.  More  details  are  given  by  Rao-Mitra  [7]  and 
Bjerhammar  [8],  among  others. 

Next  we  limit  the  maximum  step  size  for  every  component  of  Au  : 

I  Aui  |  <  dmax  (35) 

Finally  we  reverse  the  resulting  step  for  every  component,  when  it  has  the 
same  direction  as  the  gradient: 


AUi  ♦  & j  AU j 

«i  =  p(gi  Au-j  <  0)  -  p(gi  Au-j  >  0) 


(36) 


Combining  the  individual  steps,  the  following  algorithm  is  proposed: 

1.  Calculate  d  =  H*  g,  (37a) 

where  H~  is  the  generalized  Inverse  of  the  Hessian  H 

2.  Clip  the  vector  d  so  that 

d2  <  d2  (37b) 

max 

3.  Reverse  the  sign  of  d}  if 

9i  di  >  0  (37c) 

4.  Equate 

AUi  =  d1  (37d) 

The  actual  value  of  the  maximum  step  size  depends  on  the  number  of 
hydrophones.  As  a  rule  of  thumb,  one  quarter  of  the  classical  resolution 
limit  has  been  found  to  be  appropriate  for  d^x. 
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6  PERFORMANCE  ANALYSIS 

It  is  known  from  statistics  that  a  good  measure  for  the  performance  of  a 
signal -processing  algorithm  is  given  by  the  Cramer-Rao  lower  bound  (CRLB). 

Defining  z  as  the  parameter  estimation  covariance  matrix 


Z  =  E((e(y)-a)<e(y)-a)*} 
and  J  as  the  Fisher  information  matrix 


J  =  -  E{va  log  p(y  )  a)} 


where  E(...}  denotes  the  expectation  operator 

a  is  the  vector  to  be  estimated 
2 

v„  is  the  matrix  of  second  derivatives 


it  can  be  shown  [6]  that  for  an  unbiased  estimator,  the  following 
Cramer-Rao  bound  holds: 

Z  >  J-1  .  (40) 


For  the  general  case  of  multiple  sources  the  Inverse  of  the  Fisher  infor¬ 
mation  matrix  is  estimated  numerically. 


In  the  single  source  case,  however,  the  performance  bound  can  be  given  in 
closed  form  [1]. 

Let  a  be  defined  by 


a  =  (S,y  ,v) 


where 


a  = :  s  e 


1*i£(n+Yt)(l-Zu/N) 
*n,t  =:  e  x 


s  is  the  signal  (amplitude)  level 
v  Is  the  signal  (amplitude)  phase 

u  is  the  beam  number  connected  with  the  source  bearing  angle  w  via 
cos  w  *  (l-2w/N) 

N  is  the  total  number  of  beams. 
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then  from  Eq.  44 


i  2  1 

var  {b}  >  o  •  — — 
n  1280 


var  |uf  >  •  — — 

s2  4.108 


var  (v}  >  .  - - - 

s2  335.238 


(45) 


The  variance  of  the  signal  phase  will  be  ignored  in  the  detailed  analysis. 

The  performance  of  the  multiple  beamformer  is  presented  in  two  steps. 
First,  the  single  source  case  is  treated.  Here  the  standard  deviation  of 
the  bearing  and  amplitude  estimation  is  tabulated  as  a  function  of  signal 
level  to  noise  ratio  (SLNR)  and  signal  variance  to  noise  ratio  (SVNR)  where 

SLNR  =  20  log  (s/on) 

SVNR  =  20  log  (oa/an) 


The  purpose  of  this  part  of  the  analysis  is  to  obtain  an  indication  of  the 
detection  capabilities  of  the  multiple  beamformer. 

Second,  the  case  of  two  closely  spaced  sources  with  equal  amplitude  levels 
is  used  to  study  the  resolution  capabilities  of  the  processor.  In  this 
case  the  signal  level  to  noise  ratio  (SLNR)  and  the  signal  variance  to 
noise  ratio  (SVNR)  are  kept  fixed. 

Tables  1  and  2  give  the  standard  deviation  of  the  bearing  and  amplitude 
estimation  together  with  the  calculated  Cramer-Rao  Lower  Bound  CRLB.  In 
these  tables  the  values  for  the  SVNR  are  selected  as  follows: 

i)  SVNR  =  -  oo  perfect  correlation  in  time 
ii  SVNR  =  -  10  dB 
iii)  SVNR  =  -  5  dB 


The  reasons  for  this  choice  are  that 

-  in  case  (1)  the  simulated  data  fit  the  model  assumptions 

-  in  case  (ii)  the  Incoherent  part  the  signal  alone  Is  just  detec¬ 
table  with  a  conventional  beamformer  having  M  hydrophones  (for 
M  =  32  we  have  10  log  M  -  15  dB) 
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TABLE  1 


STANDARD  DEVIATION  OF  THE  BEARING  ESTIMATION 


SLNR 

CRLB 

SVNR 

_  ao 

-10 

-5 

10 

0.156 

0.181 

0.195 

0.197 

5 

0.277 

0.306 

0.313 

0.334 

0 

0.493 

0.466 

0.550 

0.521 

-5 

1.140 

0.818 

0.758 

0.845 

-10 

2.028 

1.586 

1.511 

1.395 

-15 

3.606 

2.899 

2.664 

3.022 

-20 

6.412 

6.122 

6.769 

7.515 

-25 

11.402 

11.456 

9.843 

9.369 

-30 

20.277 

17.096 

13.528 

14.789 

95%  confidence 


dB 


TABLE  2 

STANDARD  DEVIATION  OF  THE  AMPLITUDE  ESTIMATION 


SLNR 

CRLB 

SVNR 

•  oo 

-10 

-5 

10 

0.0280 

0.0244 

0.0343 

0.0533 

5 

0.0280 

0.0281 

0.0268 

0.0595 

0 

0.0280 

0.0278 

0.0335 

0.0620 

-5 

0.0280 

0.0307 

0.0312 

0.0551 

-10 

0.0280 

0.0273 

0.0305 

0.0615 

-15 

0.0280 

0.0311 

0.0286 

0.0500 

-20 

0.0280 

0.0258 

0.0352 

0.0523 

-25 

0.0280 

0.0224 

0.0285 

0.0386 

-30 

0.0280 

0.0231 

0.0261 

0.0403 

95% 

confidence 

(3:!)  - 
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-  in  case  (iii)  the  signal  can  easily  be  detected  with  a  conven¬ 
tional  beamformer  even  when  the  mean  amplitude  vanishes. 

Considering  the  bearing  estimation  (Table  1)  one  sees  that  for  all  SLNR  and 
for  all  three  SVNR  values  the  standard  deviation  follows  quite  well  the 
CRLB. 

Table  2  presents  the  performance  of  the  amplitude  estimator.  Here  also  the 
standard  deviation  is  found  to  be  very  close  to  the  CRLB,  at  least  for 
small  SVNR  values.  For  a  signal  variance  to  noise  ratio  of  -5  dB  the  per¬ 
formance  in  the  amplitude  estimation  is  4  to  7  dB  worse  than  the  CRLB. 

This  means  that  a  SVNR  of  -5  dB  is  sufficient  to  violate  the  coherence 
assumption  within  the  model  and  therefore  to  degrade  the  accuracy  of  the 
amplitude  estimation. 

Next  we  consider  the  resolution  capabilities  of  the  multiple  beamformer. 
Figures  2  and  3  plot  the  standard  deviation  of  the  bearing  and  amplitude 
estimate  as  a  function  of  the  separation  of  the  two  signals.  For  this  part 
of  the  analysis  a  signal  level  to  noise  ratio  (SLNR)  =  0  dB  and  a  signal 
variance  to  noise  ratio  (SVNR)  =  -10  dB  were  selected. 

The  two  figures  show  nearly  optimal  performance  for  the  multiple  beamformer 
for  separations  of  more  than  0.5  separation  units. 


X  N 


where 

d  =  Sensor  spacing 
X  =  Signal  wavelength 
M  =  Number  of  sensors 
N  =  Total  number  of  beams 


The  direct  comparison  with  the  calculated  CRLB  for  two  sources  shows  a 
remarkable  discrepancy  around  1  separation  unit  (half  the  width  of  the  main 
lobe),  where  the  CRLB  is  significantly  worse  than  the  simulated  perfor¬ 
mance.  These  figures  suggest  that  at  these  close  separations  the 
Cramer-Rao  lower  bound  cannot  be  used  in  the  simplified  form  of  Eq.  40. 

At  very  small  separations  the  simulation  shows  a  degraded  performance  with 
respect  to  the  CRLB.  This  is  because  the  appearance  of  outlayers  increases 
the  variance  drastically. 
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FIG.  2  PERFORMANCE  ANALYSIS:  BEARING  ESTIMATION . 


FIG.  3  PERFORMANCE  ANALYSIS:  AMPLITUDE  ESTIMATION . 
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Then  the  Fisher  information  matrix  is  given  by 


/' 


0  J  =  I 

n  ,t 


(42) 


It  is  easy  to  see  that  the  performance  bounds  can  be  written  as  follows: 


var  { s }  >  <j 


var  m  i_y _ n _ 

‘  1  s2  U  2d)  (M2 -1  )+Y2 ( T2 -1 ) 


(43) 


var  {vj  >  2_  2  (m-1)(2m-i)+3y(M-i)(t-i)V(t-i)(2T-i) 
S2  (M2-l)+y2(T2-I) 


Next  the  error  a  has  to  be  estimated.  From  Eq.  20  we  find  that: 

2 

o2  =  — 


(44) 


MT 

2 

where  o  is  the  noise  variance  as  defined  in  Eq.  4. 
n 

For  the  present  analysis  the  following  values  are  selected 


M  =  32 
T  =  40 
N  *  1024 
2d  =  A 
v  =  0 
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7  PRE-PROCESSOR 

During  the  derivation  of  the  multiple  beamformer  algorithm  it  was  always 
assumed  that  the  number  of  sources  and  their  rough  bearings  are  known.  In 
the  following  we  will  drop  this  unjustified  assumption  and  sketch  some 
possibilities  by  which  we  can  obtain  this  information. 

Let  us  first  deal  with  the  question  of  the  model  order  or  number  of  sour¬ 
ces.  Here  the  answer  seems  to  be  very  simple.  Especially  after  the  recent 
work  of  Bienvenu  and  Kopp  [9,  10]  the  favourite  should  be  the  eigenvalue 
detector. 

This  method  uses  the  eigenvalue  A,  i  =  1 . M  of  the  estimated 

correlation  matrix 


Rn,m 


\ 

T 


T-l 


l 

t=0 


★ 

yn,t  ym,t 


(46) 


to  find  the  model  order  K.  In  detail,  the  number  of  sources  is  estimated 
by  the  value  of  K,  which  minimizes  the  following  quantity 

HDL(K)  =  T  F(K)  +  \  K  (2M-K)  log  T 

with 

(M~K  M  K 

TT7  i  M)  -f  log  Ai  (47) 

i=l  /  i=l 

where  the  Minimum  Description  Length  (MDL)  after  Rissanen  is  used  Instead 
of  the  Akaike  Criterion  [10]. 

Having  now  estimated  the  model  order  we  have  to  find  the  coarse  bearing  of 
the  K  sources.  For  this  the  use  of  some  sort  of  high-resolution  technique 
is  needed  to  ensure  that  every  source  has  an  appropriate  bearing  value  to 
start  the  fine  bearing  iteration. 

Using  the  ideas  behind  the  Pisarenko  Harmonic  Decomposition  and  the  Prony 
method,  two  candidates  are  derived  [11]. 

The  basic  point  is  to  realize  that  the  signal  model  given  in  Eq.  2. 


K  nnt 

2n,t  *  I  ak  ♦k  (*8) 

k«l 
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can  be  understood  as  the  general  solution  of  the  linear  difference  equation 
K 

1  cm  zn-m,t  =  ® 
m=0 

where 

c0  =  1 
K  <  n  <  M-l 

or  in  Matrix  notation 

Zt  c  =  0  (50) 

where 


Expressing  now  the  model  in  terms  of  the  measurement  data  and  error 

Zt  =  Yt  -  Wt  (51) 

where 

Yt  is  the  measurement  data  matrix 
Wt  is  the  measurement  error  matrix 

Eq.  50  becomes  an  Autoregressive-Moving  Average  (ARMA)  model  formulation 
Yt  c  *  Wt  c  .  (52) 


Starting  with  the  ARMA  model,  two  methods  of  solving  for  the  coefficient 
vector  c  are  presented. 
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The  first  method,  which  is  connected  with  Pisarenko  Harmonic  Decomposition, 
assumes  the  coefficient  vector  c  to  be  constant  in  time  so  that  we  can  get 
an  average  ARMA  model 


(Ryy  -  o2  I)  c  =  0 
w 


(53) 


where 


Ryy  =  <Y*  Y  > 
YT  t  t 


o2  I  =  <Wt  Wt>  =  <Yt  Wt> 


<•>  is  the  time  average. 

But  this  is  an  eigenvalue  equation,  where  we  find  the  vector  c  as  the 

eigenvector  to  the  eigenvalue  o*,  which  turns  out  to  be  the  smallest 
eigenvalue  of  the  R  matrix. 

The  second  approach  takes  Eq.  52  and  replaces  the  MA  term  by  a  single  noise 
vector 


ut  =  Wt  c  .  (54) 

This  means  that  the  ARMA  model  is  modified  to  an  autoregressive  ( AR )  model 
and  the  usual  technique  to  minimize  the  prediction  error  can  be  used  to 
obtain  the  vector  c 

<  |  wt  |  Z>  ♦  min  (55) 

yields  to 

c*  Ryy  c  ♦  min  (56) 

with  the  constraint 
Cq  =  1  • 


This  approach  corresponds  to  the  Prony  method 
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The  next  step  is  to  solve  the  characteristic  polynomial  of  the  linear  dif¬ 
ference  Eq.  49 


K 

v  _  . K— m  ~ 

l  cm  ♦  =  0 

m=0 


(57) 


to  find  a  basic  set  of  solutions  {»pk{  k  =  1 . K,  which  we  need  to 

build  our  source  model. 

This  basic  set  of  solutions  is  now  used  for  the  coarse  bearing  estimates 
necessary  for  the  fine  bearing  iterations. 

However,  before  we  complete  this  section  we  suggest  a  simple  improvement  of 
these  methods. 

Spurious  harmonic  solutions  of  Eq.  57  are  suppressed  by  replacing  the  Rvv 
matrix  in  Eqs.  53  and  56  by  a  forward/backward  averaged  matrix.  This 
method  of  additional  spatial  averaging  also  gives  increased  stability  to 
the  solution  [12]. 


8  POST-PROCESSOR 

The  post-processor  serves  two  purposes.  After  having  estimated  the  signal 
components  of  the  received  soundfield  the  post-processor  can  evaluate  the 
success  of  this  estimation  procedure.  Also  it  can  show  the  spatial  distri- 
bution  of  the  signal -free  ambient  noise  power. 


The  idea  behind  the  post-processor  is  to  subtract  the  estimated  signal  from 
the  data  and  to  apply  a  conventional  beamformer  to  the  residual. 

This  means  that  when  we  subtract  the  signal  from  the  data 


wn,t  =  >n,t 


r  luk(nnt) 

l  ak  e 
k  =  l 


(58) 


and  then  estimate  the  signal-free  covariance  matrix 
T-l 


^n.m  =  j  wn,t  wm,t  =  Rn,m 

•j,"  ■  w  • . 


K  *  -iu^rn  i  T=!  -iuty t 

ak  e  ~  l  yn,t  e  + 

T  t=0 


(59) 
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we  will  get  with 


T 


T-l 

r 

t=o 


zm,t) 


iu^yt 

e 


=  0 


(60) 


the  relation 

0  =  R  -  S 


(61) 


where 


Jn 


1  T‘1  * 

~  I  zn,t  zm,t 
T  t=0 


so  that  the  ambient  noise  power  estimate  Pan  can  be  written  as 


Pan  =  v*  Q  v  =  v*  R  v  -  v*  S  v  .  (62) 

where  v  is  the  usual  plane  wave  steering  vector. 

The  first  term  in  Eq.  62  is  thereby  nothing  else  than  the  conventional 
beamformer  applied  to  the  estimated  measurement  cross-correlation  matrix; 
the  second  term  is  simply  the  signal-alone  beam  pattern. 


SUMMARY 

A  multiple  beamformer  has  been  presented  and  analyzed.  The  performance  in 
amplitude  and  bearing  estimation  is  found  to  follow  quite  well  the 
Cramer-Rao  Lower  Bound  for  varying  conditions.  The  increased  detection  and 
resolution  capabilities  are  a  consequence  of  defining  a  coherent  processor 
in  space  and  time.  However,  it  has  also  been  shown  that  small  Incoherent 
components  in  time  do  not  degrade  the  optimal  performance. 
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